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ABSTRACT 

The non-selective voltage activated cation channel from the human red cells, which is activated at depolarizing 
potentials, has been shown to exhibit counter-clockwise gating hysteresis. We have analyzed the phenomenon 
with the simplest possible phenomenological models by assuming 2x2 discrete states, i.e. two normal open/closed 
states with two different states of "gate tension." Rates of transitions between the two branches of the hysteresis 
curve have been modeled with single-barrier kinetics by introducing a real-valued "reaction coordinate" param- 
eterizing the protein's conformational change. When described in terms of the effective potential with cyclic 
variations of the control parameter (an activating voltage), this model exhibits typical "resonant effects": syn- 
chronization, resonant activation and stochastic resonance. Occurrence of the phenomena is investigated by 
running the stochastic dynamics of the model and analyzing statistical properties of gating trajectories. 

Keywords: Hysteresis, stochastic resonance, channel gating, synchronization, resonant activation, signal to 
noise ratio, spectral amplification, power spectrum 

1. INTRODUCTION 

Ion channels are proteins that allow selective conduction of ions across a biological membrane. The transitions 
between an open (conducting) and a closed (nonconducting) conformation of proteins forming the channel can 
be registered by measuring movement of ions across the cell membrane. The switching of ion channels between 
different states is brought about 1,2 by energy derived from thermal fluctuations, the electrical field and (or) 
binding the ligands. In particular, for voltage-gated ion channels discussed in this paper, kinetics between 
conducting and non-conducting states is controlled by the trans-membrane potential. The sequence of open and 
closed states can be measured in individual channels by the patch clamp technique 2 to record the ionic current 
flowing through the single channel. The resulting time series recording is a starting point for kinetic analysis 
and modeling. A common assumption for most of the models used in the study of single channel gating is that 
the protein can dwell only in a few discrete states and that the switching between these states is considered 
random with the transition probabilities depending only on the present state, and not on the previous history of 
the channel. 3 However, existing electrophysiological literature 4-9 points out existence of some degree of memory 
in the function of ion channels. The identification of short- and long-term memory in single channel recordings 
could be then an essential tool to discriminate between various models of channel kinetics and interpretation of 
the data. 

In relation to the channel memory a set of intriguing data have been reported by Kaestner et al. who have 
studied 6 voltage-activated cation channel in the human red blood cell (RBC) membrane, demonstrating that 
the degree of activation at a given membrane potential depended strongly on the prehistory of the channel. The 
gating exhibited hysteretic-like behavior with the quasi steady-state deactivation and activation curves displaced 
by about 25 mV. The effect has been explained 10 by evaluating a kinetic model for voltage-activated channel 
gating that leads to a delayed response to a cycling external field in a way similar to the change of magnetization 
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in ferromagnetic materials under the influence of an external magnetic field. Alike memory-phenomena have 
been also discovered in the ligand-gated NMDA receptor, 11 in voltage-gated HERG-like channels 12 and in the 
gating kinetics of the vanilloid receptor 13 which is a thermo-sensitive non-selective cation channel activated by 
noxious heat. In order to explain hysteretic behavior observed in channel response to external stimuli some of 
the authors 12 suggested the existence of an ultra-slow inactivation process that effectively introduces a lag in the 
feedback between voltage and gating. Alternatively, 13 a periodic sensitization of the channel on repeated heat 
stimulation has been assumed. 

Viewed as a nonlinear and time-delayed response of a system to the cyclical variation of a control parameter, 
hysteresis is a fairly common phenomenon observed in bistable dynamical systems. 14 ' 15 In particular, it has 
been argued 16-18 that the hysteresis loop area can be a useful quantity for identifying critical values of field 
parameters responsible for the onset of resonance-like response of the system and therefore helpful in designing 
field parameters for an optimal control. In the series of articles, Mahato, Shenoy and Jayannavar 15, 17 performed 
studies aimed to understand connection between the stochastic synchronization and stochastic resonance (SR) in 
a bistable system exhibiting hysteretic response to the external saw-tooth-type driving. They have shown that 
the hysteresis loop defined in such a system is a good measure of synchronization of passages from one well to 
the other under the influence of a Gaussian white noise. Moreover, at a particular value of the noise intensity, 
the maximum of the signal-to-noise ratio (SNR) is achieved thus demonstrating appearance of the stochastic 
resonance, that in turn, can be correlated with the hysteresis loop area. These findings suggest that hysteretic 
voltage-dependent gating may be closely related to documented observation of stochastic resonance in biological 
channels. 19 

Since the discovery of the SR phenomenon more than 20 years ago, the idea of stochastic resonance has gained 
a lot of interest in diverse disciplines of science where its presence has been demonstrated repeatedly 20 in versatile 
applications. At the subcellular level, existence of the SR has been also proven in the ion channels 19 where the 
sensory transduction of the signal (ionic current) can be optimized by the presence of intrinsic ambient noise. The 
result has provoked a question whether the SR is a collective effect occurring in some finite assemblies of channels 
or, perhaps could be also realized in a single channel activity. Further theoretical studies have revealed 21 that in 
a single Shaker potassium channel the SR phenomenon can occur if the parameters that describe its functionality 
are set within the regime where the channel is predominantly dwelled in the closed state. 

In this contribution we discuss occurrence of resonant effects in a model double-well system that has been 
used to describe voltage-dependent gating kinetics in the RBC cation channel. 10 Section II presents basics of the 
reaction pathway modeling for that case. The interplay of voltage driving and thermal fluctuations are studied in 
an overdamped limit and lead to hysteresis, stochastic resonance and resonant activation. All those phenomena 
can be registered in the analysis of the channel gating kinetics, as discussed in Section III. Section IV is devoted 
to a standard two-state approximation within the Kramers scenario. Summary and closing remarks are presented 
in Section V. 



2. BISTABILITY AND REACTION PATHWAY 

The mechanism of voltage-gated ion channels was introduced in an early work of Hodgkin and Huxley who 
suggested that Boltzmann statistics describes their probability of being open or closed, including the voltage- 
dependence of this probability. 1 ' 2 For the purpose of modeling, we assume the ion channel in question has only 
two states, open and closed. We associate an internal energy E opcn , -E c ioscd with each of the two states, and 
assume that Boltzmann statistics describes the probability for a channel to be open or closed, with transitions 
between the two states being thermally activated barrier crossings in an energy landscape that we do not detail. 
Then the probability that a channel is open (p op en) or closed (p c iosod) is 

g-^open, (closed) i 
Popon, (closed) = e _ Bopon + e -B cloaod = 1 + e ±AE ' W 

where all energies are given in units of k^T. In the case of voltage-gated channels, the transition between the 
open and the closed state possibly involves movement of charges in the trans-membrane field, reorientation of 



local polarization and conformational variation of channel protein that results in the change of energy AE of the 
system. With AE as above, Eq. (1) may be written 10 

Po^(V)= i + e _ 1 a{v _- ) =f(V), (2) 

where aV replaces (-Enclosed— E open )/ltBT, V is the voltage at which p opcn = p c iosod = 1/2 and ceV parameterizes all 
other contributions to AE. Equation (2) may be derived by assuming Kramer's scenario for the gating kinetics. 22 
This approach offers a description in which the configurational change of the channel protein is modeled by a 
multi-component "reaction coordinate" x. Since the potential energy of the entire system (the gate, channel 
protein and membrane environment) is a function of hundreds of relevant coordinates of the system, x can be 
understood as an effective variable similar to its definition in the Marcus theory of charge-transfer reactions. 23 
The open and closed states of a channel can then be described as x w 1 and x « — 1, respectively, with a potential 
(activation) barrier separating these two regions on the x-axis, and x being confined to the regions near —1 and 
1 by the potential taking large values elsewhere. Opening and closing of a channel is then modeled by a process 
similar to an activated uni-molecular chemical reaction pictured by a kinetic scheme 

feo(t) 

open f± closed. (3) 

fco(t) 



In order to capture the configurational change of the protein quantitatively, we introduce a parameterization for 
the reaction coordinate x for this change by assuming that an equilibrium probability for the system to be in a 
conformational state x is given by 

(4) 



-W e tt(x;V) 



CT=±1 



where a describes a binary variable a — ±1 discriminating between open (tr = 1) and closed (<r = —1) states of 
the channel. A corresponding effective potential (cf. Fig. 1) takes on the form 10 



W eS (x;V)=W(x)-log 



cosh ( i(A£ gato + AE int x] 



where W{x) describes the protein's internal energy as a function of x 

(\ ° 
W(x) = AE con{ -x(3 - x 2 ) + 



16x n 



-(l-x'f 



(5) 



(6) 



with A_E con f being a measure of the difference in the internal energy associated with the change of protein's 
configuration favoring open (closed) states of the gate. The parameter x max has been chosen as the value of 
the reaction coordinate x at which the local maximum, the top of the barrier between x = ±1 states is located; 
AE gatc — — a(V — V) is contribution to the internal energy discriminating between the open (closed) states of 
the gate and A_E; n t = —aAV describes interaction energy between the gate and surrounding protein. In the 
above formula AV stands for the maximum intensity of driving external voltage V and V is an average value of 
V for which stationary probability p pen achieves value 1/2. 

We now additionally assume that the membrane potential felt by a channel is not exactly constant, but 
changes, albeit slowly so compared to the rates of the gate opening and closing 

V(t) = V- Ay cos fit. (7) 

For V(t) — V large positive, stability analysis shows that conformational state x = 1 is favored over x = — 1 with 
an energy difference aAy = |A£i nt |. Similarly, for V(t) — V large negative 



W cS {x; V) = W(x) + ^a(V — V + AVx) + log 2 



(8) 




Figure 1. The effective potential for conformational changes, W e s(x; V), for V = V, V = V + 3AV and V = V — 3AV. 
Parameters of the potential (in units of k B T) AE cont = 0.6, A£ gatc = 9.72, AE int = 2.16, x max = 0.018, AV = 12 and 
V = 66. The set of parameters has been chosen to model hysteretic behavior observed in opening probabilities in the 
RBC channel. 6 

and the barrier at ir max becomes lowered by \ct\V — V\ = ^|A.Eg a te(V)| favoring state x = — 1 over x = 1. Note 
that for large differences V(t) — V, addition of time-dependent driving does not change characteristic relaxation 
frequencies of the system within the wells of the potential W e g(x). Direct evaluation of = W" s (x niax ) 

and Wj.! = W" s (x = ±1) yields w^ max « -3A£' c /4, lo^ w 3AS c /(2x max ) for values of x max close to 0. We 
further assume that fluctuations of the collective variable x quickly equilibrate, i.e. their dynamics can be well 
approximated by an overdamped Langcvin equation 

dx _ dW cS (x;V) 
dt dx 

with representing effects of thermal noise modeled by uncorrelated Gaussian white noise, i.e (£(t)£(t')) = 
<r 2 S(t — t 1 ). For convenience, we set the friction coefficient 7=1 that rescales the time t in the above equation. 
Fig. 2 presents samples of stochastic trajectories X(t) simulated from the above Langevin equation at constant 
intensity a 2 = 8 of the additive noise (mimicking conditions of a constant temperature of the thermal medium) . 
They model records of current in a biological channel switching between conducting and non-conducting confor- 
mations. As expected from the former analysis, at large positive values of the gating voltage V, the channel is 
predominantly open with current fluctuating due to the local conformational variations around the open state 
(x = 1). At V = V none of the states is favored, whereas for large negative V a preferential stay in a closed state 
is visible. Increasing the noise intensity at a constant value of the voltage cycling frequency causes a trajectory 
to fluctuate wildly between both wells of the effective potential shortening the time the process spends in the 
vicinity of the metastable points x = ±1 and eventually, destroying the pattern of the signal switching between 
"open" and "closed" states. 

From the residence time frequency histograms for the open and closed states, we have replotted probability 
of the channel being open as a function of the driving voltage. Fig. 3 presents estimated opening probabilities 
for various intensities of thermal noise. 

By increasing the activating voltage from V = VtoV = V + 3AV and reversing the direction of voltage 
changes to V = V — 3AV, an apparent hysteretic shape of the p pen(V) curves is observed (cf. Fig. 3). The 
hysteresis changes its shape and area with increasing intensity of the thermal noise and disappears for large 
values of a 2 . Through cycling V, the height of the potential barrier dividing wells of open and closed states is 



(9) 



V=V+3AV 
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Figure 2. Sample realizations of the stochastic process X(t), for the two state and the continuous model, derived from 
numerical simulations of Eq. (9). Time step of simulations has been set up to 1CP 3 with the frequency of the cycling 
voltage 57 = 1. Trajectories have been run for the duration time T « 10-7T. Initial value of X has been sampled from the 
uniform distribution on the interval [—1.5, 1.5]. 




Figure 3. The open state probability p pen as a function of the membrane potential V. Probability has been estimated 
from the frequency histograms for 3 x 10 3 trajectories mimicking passages between different conductance levels. Intensity 
of the thermal noise a 2 = 6 (left panel) and a 2 = 8 (right panel); time step of the simulation dt = 10 -3 . 
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Figure 4. Power spectra of the channel gating process for the two-state (left panel) and continuous model (right panel), 
as simulated from Eq. (9) at the intensity of the additive thermal noise a 2 = 0.2 and Q = 1. 

now itself a time-dependent variable. From other contexts, 19 ' 20, 24, 25 it is known that this presence of periodic 
driving force and additional noise (fluctuations) may lead to several nontrivial resonant effects which we are now 
going to address in the forthcoming section. 



3. DYNAMICAL HYSTERESIS, STOCHASTIC RESONANCE AND RESONANT 
ACTIVATION IN CHANNEL GATING KINETICS 

An essential property of the stochastic process {X(t)} is interdependence of the random variable X measured by 
the covariance function C(t 7 r) = (X(t + r)X{t)) which through the Wiener-Khinchin theorem is related to the 
power spectrum of the process {X(t)}. In the framework of Fourier transformation, properties of the covariance 
function may also be read off from its Fourier transform for a discretized version of trajectories X(t). Figs. 4 
and 5 present plots of the power spectral density as a function of frequency from a simulation of Eq. (9). The 
power spectra have been obtained from discrete periodograms for two different values of the intensity of the 
additive noise At low noise intensity a 2 = 0.2, the power spectrum exhibits sharp maxima at multiples 

of the unit frequency of the driving voltage V (cf. right panel of Fig. 4) . The existence of spikes in the power 
spectrum is a general feature for periodically driven systems 26 and reflects (asymptotically) periodic character 
of the probability distribution p(x, t). This feature is absent in the power spectrum of the two-state model (see 
left panel of Fig. 4). In this case any trajectory reaching the top of the potential barrier at x max is assumed to 
switch the system from one state to another, so that the outgoing "signal" X(t) is filtered and stored as a binary 
series of ±1 values. On increasing the noise intensity (or, cquivalently temperature of the channel environment) 
to (j 2 == 6, the structure of multiple maxima washes out leaving the dominant peak at u> « fl = 1 observable for 
both - the two state approximation and the continuous model (cf. Fig. 5). For large values of uj, power spectra 
scale to the uu~ 2 form reflecting exponentially decaying memory of the gating process {X(t)}. 

The spectral amplification (an integrated power in the time-averaged power spectral density at w « = 1) 
has been calculated (cf. Fig. 6) according to the interpretation due to Jung and Hanggi. 26 The mean value of the 
{X(t)} process has been obtained by averaging over the ensemble of various realizations of the additive noise. 
For long times t^ooan asymptotic trajectory looses memory of the initial conditions and becomes a periodic 
function of time. Its square amplitude is then proportional to the spectral amplification expressed as the ratio 
of integrated power stored in the delta-like peaks of the power spectra at ±fl to a total power of the modulation 
signal in V. 




Figure 6. Spectral amplification, in arbitrary units for the continuous model, of the signal (registered as a trajectory 
passing from the left to the right well of the potential W e g (x; V)) for f2 = 1 at various noise intensities a 2 . 
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Figure 7. The signal-to-noise ratio for f2 = 1 for the two state model (left panel) and for the continuous model (right 
panel), as a function of the noise strength a 2 . The decaying part of the curve, in the right panel, represents relaxation 
processes in the potential well, details in the text. 



Spectral amplification displays a typical bell-shaped stochastic resonance behavior as a function of increasing 
noise intensity with a characteristic peak observed for a 2 w 6 for both models under study. Similarly, evaluation 
of the signal-to-noise (SNR) ratio for the potential under study produces curve (Fig. 7) that exhibits non- 
monotonic resonant behavior for a 2 « 6 and characteristic divergence for the noise strength a 2 — ► 0. For small 
a 2 (cr 2 < 2) barrier crossing events are rare events, and the decaying part of the SNR curve (cf. right panel of 
Fig. 7) corresponds to relaxation processes within the potential well. This feature is absent for the two state 
model approximation. In accordance with the former studies of Mahato and Jayannavar, 15 ' 17 the maximum 
area of the hysteresis curve is observed for the noise intensity at which the stochastic resonance is registered. 

The same system may also exhibit a different resonant behavior - so called resonant activation (RA) 25,27 
expressed by the shortest time of passage over a fluctuating potential barrier in the presence of an additive noise. 
As expected based on previous theoretical studies, 25 ' 27 ' 28 the RA occurs typically in the bistable system under 
broad circumstances of varying shape of the potential barrier and qualities of fluctuations. In the case of the 
modeled channel system, the resonant activation is brought up by an interplay of cyclic (deterministic) variations 
of the barrier and an additive thermal noise (Fig. 8). 

In order to examine the statistics of the first passage times FPT for various shapes of the potential W e s(x; V) 
introduced by the cycling voltage V(t) own program performing Monte Carlo simulations has been used. The 
trajectories Eq. (9) have been run starting from the left well until the top of the barrier has been hit, where 
by the assumption, the trajectory has been absorbed. A position of the top of the barrier x max depends on the 
driving voltage (see Fig. 9) and for a given set of the parameters smoothly changes between x max « 0.06 and 
Xmnx ~ —0.03. Within the error limits of the simulations the resonant activation becomes visible at the noise 
intensity a 2 — 2 and the effect vanishes with increasing values of a 2 . 



4. TWO STATE APPROXIMATION 

In order to discuss the configurational change of the protein quantitatively we have introduced a parameterization 
of the reaction pathway for this change. We now explore the consequences of the specific form of W e ff(x; V), 
approximating the thermal barrier crossing with the simple Kramers scenario. In doing that, our description 
captures, albeit approximatively only, the consequences of a much larger set of potentials which are generic for 
studying bistability and SR phenomena. 
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Figure 8. The MFPT for trajectories starting from the left minimum of the potential W c g(x; V) as a function of driving 
voltage frequency f2 at the noise strength a 2 = 2. Numerical results were obtained by use of Monte Carlo simulation of 
Eq. (9) with time step dt = 1CP 5 and averaged over N = 10 3 realizations. Error bars represent deviation from the mean. 
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Figure 9. Location of the barrier top for periodically perturbed potential as a function of the driving voltage. Parameters 
AS con f, ASi„t, A£; ga t , AV and V as in Fig. 1. 



In the overdamped Kramer's scenario, 29 after assuming that the thermal energy k^T is much smaller than 
the respective barrier height A14 71 * = W c g (x max ; V) — W c e (±1; V), the reaction coordinate undergoes a creeping 
motion and an adiabatic elimination of the velocity x leads to the Langevin equation (9). The time evolution of 
the corresponding probability p(x, t) is governed by the Smoluchowski equation. The spatial-diffusion rate in the 
potential W e s(x;V) can be then evaluated by a steepest-descent approximation from the quadrature formulas 
for the mean first passage time (MFPT) 

k = k-i^+i(V(t)) 

= ^_ e -(Weff(x„ax;V(t))-W eff (-l;y(t))) ; 

= v+e -(We a (x^;V(t))-w aa (+i,v(t))) ^ (20) 

where the subdominant factor u± = (lu ±\u> Xmax ) /2n contains product of the characteristic angular frequency 
inside the metastable minimum (x = ±1) and an angular frequency at the top of the potential barrier. 

For a sufficiently high barrier separating open and closed states, the probability of x taking values in between 
vanishes and we can work with the integrated probability of x being near —1 or 1, n c i ose d(t) = f*™* p(x, t)dx 
and n opcn (t) — +s p(x,t)dt where {—6,5} describes an interval in the vicinity of the metastable x max . The 
integrated probabilities satisfy n c \ osc( i(t) + n opcn (t) — 1 and the linear kinetic equation 

— = fco'T-closod(i) — fcc^opon(^)- (H) 

Under stationary conditions the ratio of the reaction rates can be related to the difference in the energy by the 
equilibrium constant 

K _ K _ "closed _ ^+ e -[Weff(-l;V(t))-Wef f (l;y(t))]_ ( 12 ) 
^open V— 

Rearrangement of this equation gives a formula for the probability of the channel being open Eq. (2) 

Popcn = i + k c /k Q - (13) 

Thus AWeff is related to the ratio of kinetic rates for the "forward" and "backward" reactions. Experimental 
estimates of the opening and closing rates can be found from the voltage-dependent dwell times determined 
from experimental recordings as inverse of averages i.e. k = (T c ) _1 , k c = (T ) _1 where (T C;0 ) are average dwell 
times in the closed (open) states. The gating dynamics described by Eq. (11) can be otherwise considered as a 
two-state Markovian process 30-32 with a time-dependent transition probabilities between the states. A two state 
model of that type is known to describe a stochastic resonance 20,33 and with both periodically and randomly 
modulated rates has been investigated.. 32 

The general solution of Eq. (11) can be readily found 

^opcn(^) — g I c( j ) o( ^Tl O p en (0) 

+ C f dsk (s)e- [K ^ s)+K ^ t - s) \ (14) 
Jo 

where t 

Ko,c(*,s) = I dt'k , c {t') (15) 

J s 

and 

k .c = K 0tC e~ a °' cV [1 + a 0:C cost + A 0iC r/(i)], (16) 



where rj(t) stands for other sources of (nonthermal) noises of intensity A 0jC that can perturb the transition rates. 
Note, that the Kramers rate formula is derived under the assumption that the probability density within the 



well is close to its equilibrium distribution, so that in order to use the modified Kramers rate with a periodic 
perturbation of the voltage, its frequency must be much slower than the characteristic rate for the probability 
to equilibrate within the well. Such an adiabatic approximation is valid if the second derivative of the effective 
potential W c h(x;V) at x = ±1, i.e. u>±i is much bigger than 1 (in our analysis the characteristic period of 
the cycling voltage has been set up to 1). Parametric plots of the cumulative open probability, as defined by 
^opon in function of V display then the dynamic hysteresis whose shape and loop-area depend on the ratio of 
the driving voltage frequency to the overall relaxation rate k = k Q + k c . The hysteresis disappears only in case 
of driving that has time scale significantly different from the inverse of rates k 0tC , i.e. when it is either very slow 
or very fast. In such situations, the response of the system measured in terms of the integrated probabilities 
ftopen, closed depends on the applied voltage but not on its frequencies (slow driving), or follows (fast driving) 
average, time independent rates. In consequence, hysteretic behavior may be considered a generic for systems 
exhibiting the SR phenomenon. For the periodically driven two state process, recent analysis performed by 
Talkner 34 has shown that studies of the entrance time densities in the system can detect the synchronization 
of the driving force and the driving process. However, a more detailed comparison of stochastic resonance and 
synchronization 17 ' 34 shows that the effects are closely related but not identical. 

5. SUMMARY 

There is no unique consensus in the SR literature on how to measure the SR effect. The most popular performance 
measures used in the study depend on the forcing signal and noise and can vary from system to system. The 
common approach to characterize SR is the SNR that measures how much the system output contains the input 
signal frequency 0. An adiabatic approximation 20 ' 33 can then give an explicit SNR ratio for the quadratic 
bistable potential. Alternative ways to measure SR are related to the cross-correlation measures and estimates 
of the density probability of residence time 17 and escape rates. 34 

In this work we have concentrated on a class of resonant effects that can be exhibited in the gating kinetics 
of the voltage-driven biological channel. Our analysis reproduces all typical phenomena that are characteristic 
for the two-state system driven by fluctuations and periodic forcing. Albeit the latter does not function here as 
an additive periodic force but rather modulates the height of the free energy barrier that discriminates between 
different conformational states of the gate, its presence and interference with inherent thermal fluctuations leads 
to a fine-tuning of the gating kinetics exhibited by stochastic resonance and resonant activation. 
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